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We investigate the dynamic quantum tunneling between two attractors of a mesoscopic driven Duffing oscil- 
lator. We find that, in addition to inducing remarkable quantum shift of the bifurcation point, the mesoscopic 
nature also results in a perfect linear scaling behavior for the tunneling rate with the driving distance to the 
shifted bifurcation point. 
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I. INTRODUCTION 

In the light of nanomechanics flj] and as a qubit readout 
device for superconducting qubits (i.e., the Josephson bifur- 
cation amplifier) J2-Htl, the quantum dynamical properties of 
the driven Duffing oscillator (DDO) gained renewed inter- 
est in the past years @-[I2l> since these real systems can be 
well described by a DDO under proper parameter conditions. 
While some experimental results can be understood by classi- 
cal means ILlUlql . the success of such interpretations is ba- 
sically due to the fact that the quantum levels are hardly dis- 
tinguishable from classical nonlinear excitations 11411 . and the 
data were classical in origin flfjj . Actually, possible quantum 
behaviors and the typical features of the DDO are of great in- 
terest to the community mentioned above. For instance, in 
quasiclassical limit, the quantum feature in the bistable re- 
gion is investigated against the conventional classical DDO, 
by simulating a Lindblad-type master equation and compar- 
ing the Wigner function with classical probability distribution 
in phase space Interesting result reached there reveals 
that the quantum effect makes easier the transition between 
the bistable states. Moreover, in Refs. Isl- HlTl . the switching 
rate between the bistable states near the bifurcation point, ow- 
ing to quantum and/or thermal fluctuations, were carried out 
by means of the WKB theory or semiclassical methods such 
as the mean-first-passage-time approach. 

While much more work should be done for unambiguous 
demonstrations of the quantum behaviors of the DDO, stud- 
ies of the quantum dynamics of DDO in certain fully quantum 
regimes are valuable and necessary at the early stage before 
attacking the whole difficult problem. Indeed, in a fully quan- 
tum regime which only involves a few quantum levels of the 
potential well into the dynamics, quantum behaviors of the 
DDO such as resonant tunneling and photon-assisted tunnel- 
ing are discussed in terms of amplitude and phase responses 
to the driving frequency |7|]. In the present work, we consider 
further an intermediate quantum regime, say, the quantum dy- 
namics of a mesoscopic DDO, which involves more than ten 
quantum levels in the nonlinear dynamics. Obviously, this 
is a regime between the quantum few-level and the classical 
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dense-level (or continuum) limits lfl7ll . where the quantum ef- 
fect should be crucially important, and at the same time the 
DDO's nonlinear characteristics can emerge as well. More 
specifically, we will focus on the dynamic quantum tunneling 
between the two attractors in the bistable region ifTill . in par- 
ticular accounting for a quantum shift of the bifurcation point 
associated with the mesoscopic nature, and extracting a new 
scaling exponent with the driving distance to the shifted bifur- 
cation point. 

II. MODEL AND QUALITATIVE CONSIDERATIONS 

The DDO is typically described by 

H s {t) = p 2 /2m + mQ. 2 x 2 /2 - yx A + F(t)x, (1) 

where F{t) = 2Fq cos(vf) stands for the external driving force. 
Moreover, a DDO is inevitably influenced by its surrounding 
environment, which can be quantum mechanically modeled 
by a set of harmonic oscillators, He = 2;( m ;or.^/2+pr/2m;). 
In a spirit of the Feynman-Vernon-Caldeira-Legget model 
lfl9ll . the coupling between the Duffing oscillator and the en- 
vironment is well described by Hj = -iX/^rt- After in- 
troducing the spectral density function, J(u>) = n £j A 2 5(co - 
w/)/(2m/W/), an important situation is the so-called Ohmic 
case, corresponding to J(a>) = mKQ)exp(-u)/ u) c ), where k is 
the friction coefficient, and a> c a high frequency cutoff. In the 
Ohmic case, one can successfully recover the classical fric- 
tion force lfl9tl . which plays an essential role in the classical 
dynamics of DDO 0]. In this work, we will employ a mas- 
ter equation approach to simulate the quantum dynamics of 
DDO. This treatment is equivalent to the (quantum) Langevin 
equation in which not only the friction but also the fluctu- 
ations are taken into account. For the convenience of later 
use, we also introduce here the operator b — 2, /I,/?,/ V2, with 
b, — (mjCJiXi + ipi)l ylmficoi. 

Notice that the Duffing oscillator described by Eq. (HJ has 
only finite number of bound states. In the absence of exter- 
nal driving, this becomes clear from the so-called soft non- 
linearity (y > 0) potential profile @], V(x) = ^x 2 - yx 4 , 
which defines a single well with identical barrier height Vo = 
m 2 Q. 4 /(I6y) at x = ± ^mQ. 2 1 (Ay). Nevertheless, compared to 
the hard nonlinearity (y < 0) potential well 0], which sup- 
ports infinite number of bound states, the resultant nonlinear 
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dynamics has no fundamental differences. This is because we 
are interested in a region not too far from the bottom of the 
anharmonic potential well. For the soft nonlinearity potential 
well, the chosen parameters can assure that the unbounded 
states have negligible effect on the dynamics. This also in- 
dicates that the "metastability" to be addressed in the follow- 
ing has no relevance to the "unbounded" states at very high 
excitations in the soft well. As a rough estimate, the num- 
ber of bound states is the ratio of Vo and HQ., which gives 
N = T0m = = m- Here we introduced S = mCl/h, 
and y = y/(mQ. 2 ). In our model, y = mQ 2 /24, so approxi- 
mately the number of bound state is 3S/2. In the experiment 
of Ref. JH, for instance, a rough estimate gives S 366. In 
the present work, however, we will set N ^ 12 by altering the 
circuit parameters. This defines a mesoscopic regime for the 
DDO with 18 quantum levels involved in the dynamics. Also, 
we assume an accessible temperature of 5 mK. 

Now we present a qualitative analysis to the driving dynam- 
ics of DDO, respectively, in the laboratory frame and in a ro- 
tating frame. In Fig. 1(a) we show the energy level diagram of 
the Duffing oscillator in the absence of driving. To the second- 
order perturbation of the quartic potential, the energy level 
reads E„ = [«+ 1/2- 3y(2n 2 + 2n + 1)/(4N)]MX Accordingly, 
the adjacent level spacing, AE„ = E„ —E n -\ = (1 -3yn/K)Ml, 
decreases with «. This property, together with a negative fre- 
quency detuning (e.g. 5 = 1 - v/O = 0.065 in later simu- 
lation), would result in the bistability behavior. Qualitatively 
speaking, for weak driving, the oscillator will largely remain 
in the initial ground state; for stronger driving, however, it 
will be increasingly excited to high energy states around n*, 
roughly determined by Kv = AE n *, leading to n* = N5/(3y). 

We notice that an alternative (a better) way to understand 
the quantum dynamics is in a rotating frame with the driv- 
ing frequency v, where the driving field becomes time in- 
dependent. This can be implemented by applying a rotat- 
ing transformation U{t) = exp{-;vfa'a}, where a (cr) is 
the annihilation (creation) operator of the DDO. Dropping 
fast rotating terms, i.e., under the rotating wave approxima- 
tion (RWA), we obtain H s = 6[p 2 /(2m) + (l/2)mQ 2 x 2 ] - 
6y/(4m 2 n 4 )[p 2 /(2m) + ( 1 /2)mfi 2 x 2 ] 2 + F x. In the absence of 
driving, this rotated Hamiltonian is naturally diagonal, with 
the eigenstates of harmonic oscillator iff,,, and eigenvalues 
E„ = [n + 1/2 - 3y/(2N<5)(« + 1/2) 2 ]M£2, as schematically 
shown in Fig. 1(b). Interestingly, we find here that the n*th 
level is the highest one (noting that n* labels the resonant 
level-pair in the laboratory frame). Also, the level spacing 
in the rotating frame is much smaller than its counterpart in 
the laboratory frame. Then, it is clear that the DDO's dynam- 
ics in this rotating frame is governed by the interplay of Fqx- 
induced transition and environmental dissipation. More inter- 
estingly, in the presence of Fox, we can diagonalize the trans- 
formed Hamiltonian Hs, and denote the eigenstates by iff,,. 
It is found that a one-to-one correspondence exists between 
iff,, and iff,,, which can be understood in the spirit of adiabatic 
switching. Then, we may regard i/v- and (//^-dominated states 
(wavepackets) as two attractors, which direct the evolution, 
determine the final state, and are responsible to the bistable 
behavior. Here we may mention that the classical dynamics 
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FIG. 1: (Color online) Energy diagram of the non-driven Duffing 
oscillator in (a) the laboratory frame, and (b) the rotating frame. 



of DDO can lead to possible chaotic motion, see, for instance, 
a brief discussion in Ref. 14J]. For the purpose of quantum 
detector applications, however, the chaotic motion should be 
avoided by setting a low bound of detuning J3l ■ In the present 
work, as we did previously in Ref. rf20ll . we have similarly set 
an upper bound of detuning in order to eliminate the possible 
chaotic motion. In the mesoscopic regime, the iff,,*- and \J/q- 
dominated wavepackets, which can be regarded as the coun- 
terparts of classical attractors, consist of a couple of discrete 
states. In the following, these "quantized" attractors will be 
named also small amplitude state (SAS) and large amplitude 
state (LAS). 



III. MASTER EQUATION IN ROTATING FRAME 

To account for the environmental effects on the DDO, a 
Lindblad-type master equation can be employed fiSESISl^ 
simply following the typical treatment for cavity photons (os- 
cillators) in quantum optics. In this work, however, viewing 
the "x-xf type interaction assumed by the Feynman-Vernon- 
Caldeira-Legget model which is actually beyond the rotating 
wave approximations iH HI, [1(3, [HI], we would like to use a dif- 
ferent type of master equation, say, the Born-Markov-Redfield 
(BMR) version 12 1142311 . to study the quantum dynamics of a 
weakly damped DDO as it is. Following the notation proposed 
in Ref. I22I1 . the BMR master equation for DDO's reduced 
density matrix reads 

p(t) = -ifr 1 X s p(t)-H- 2 f driZAWit^L^it^pit). 
Jo 

(2) 

Here -C s ( - ) = [H s , (•••)] is associated with the DDO 
Hamiltonian in laboratory frame, while @(t, t)(- ■ ■ ) = 
G{t, r)(- • • )G' (f, t) is the propagator in Liouvillian space. The 
interaction Liouvillian X/ is defined by £,/(■ ■ ■ ) = [Hj, (■■■)], 
where Hi = —x 2 j XjXj = -xXe describes the coupling 
of DDO to environment. The average in Eq. © means 
((••■)) - tr £ [(- ' ' )Pe\, with Pe the thermal equilibrium 
density operator of the environment. In laboratory frame, 
the driving in Hs is time-dependent, which would compli- 
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cate the dissipation terms in Eq. (O and make the numerical 
simulation difficult. To overcome this difficulty, we trans- 
form Eq. (ffj) into the rotating frame: pit) = -ihr l £.sp(t) - 
fT 2 ^dT{Zi{t)@(t,T)ti(T)g\t,T))p(t). Here, the various 
transformed quantities are defined as: pit) = U^p(t)U, 
£ s (---) = [rfHsU + ihtfU, (■■■)] = [# s , (•■■)], &(••■) = 
[U^HjU, (■■•)] = [Hi, (• • ■ )L and @(t, r) is the propagator as- 
sociated with Hs . Note that in the rotating frame the coupling 
Hamiltonian becomes time dependent, i.e., Hj = -U^xUX E = 
- y/h/2mD.(a' e ,vt + ae~' v, )X E . Inserting this result into the 
EOM of p(t) yields 

Pit) = -ih- l ls~pit) - ^{[a\ id-Is + v)a)p] 

+ [a, iCi-ls - y)a')p] + e i2v '[a\ (C(-£s ~ v)J)p] 
+ e- i2v '[a, iCi-£ s + v)a)p] + H.c.}. (3) 

In deriving this result, the well established Markov-Redfield 
approximation has been applied. Accordingly, the spec- 
tral function, C(Xs), is a Fourier transform of the environ- 
ment correlator: CiZs) = dtCit)e m ' llst , and C(t) = 
Tr E [X E it)X E iO)p E ]. 

Eq. (0 is the desired equation we obtain in the rotating 
frame. To relate it with other work, we first drop the fast oscil- 
lating terms in Eq. (f3]l under the rotating wave approximation. 
Then, we assume further a harmonic approximation: ~ 
-h5Q.a, and X.sa' « HSCla*. As a result, the well-known 
Lindblad master equation is obtained: pit) = —ifT 1 X.spit) + 
k{[1 + n(Q)]D[a]p + n(£l)D[a*]p] , where the Lindblad su- 
peroperator is defined through D[A]p = ApA' - \{A T A,p}. 
In obtaining this result, the explicit form of the spectral func- 
tion has been used, i.e., C(w) = 2[1 + n(oj)]J((o), where n(w) 
is the Bose function. We notice that in Ref. iflQIl the study 
of quantum activation was actually based on this same equa- 
tion, together also with a few techniques and approximations 
involved in the later analysis. 

IV. QUANTUM SHIFT OF THE BIFURCATION POINT 

For comparative purpose, we first outline the result from 
a classical analysis. It is well known that the classical DDO 
obeys fl]: mx + m£l 2 x + micx — 4yx 3 = -Fit). Define dimen- 
sionless variables r = Qf, w = v/Q, / = yj4y/im 3 Q. 6 )2Fo, 
Q = Q//c, and A = -2Q5; and introduce rotating transforma- 
tion x(t) = [xir)e iwT /2 + c.c.]/ V4y/(mQ 2 ). Then the slowly 
varying amplitude x(t) of the DDO in the rotating frame satis- 
fies the following EOM fl]: 2/x(t) = [iA-i)/Q+3\x\ 2 /4]x-f. 
Stationary solution of this equation yields a bistability dia- 
gram as shown in Fig. 2(a). Moreover, two critical driving 
strengths, /b(A) and f B iA), can be obtained in the limit Q » 
1: f BrB iA/A c ) = f c /i2x^ 2 ) yfl + 3x 2 ± (1 - x 2 y/ 2 \ x=Ac/A . 
Here A c = - a/3, and f c = 2 5/2 /(3 5/4 Va 1 ). Accordingly, after 
restoring dimensional units, the critical driving strengths read 
Fg g — -y/m 3 Q 6 /( 1 6y)f B B . In what follows, we will focus on 
the most important upper point Fg, which is to be re-denoted 
as F c and taken as the unit of the driving force. 



Unfortunately, in mesoscopic regime, as we will see later, 
the critical driving strength F c determined above does not 
match the result from numerical simulation. It would thus be 
desirable to develop a quantum version for the EOM of x(r), 
in order to reach consensus with the direct numerical simu- 
lation. Following Ref. l24ll . given that the reduced density 
matrix satisfies the Lindblad master equation, in Heisenberg 
picture the operator a should obey an equation of motion as 
follows 

ait) = -itT y [a,H s ] + KD[a]a = -ifr l [a,H s ] - -a, (4) 

where the ^MaZ-Lindblad superoperator is defined through 
D[A]a = A^ aA-\{A^ A,a). Here we assumed the low temper- 
ature limit «(Q) <k 1 . Moreover, below F c and starting with 
a small-amplitude state, the subsequent evolution will largely 
remain in a coherent state I25I - I27I1 . Then, in coherent state 
representation and relating the coherence number air) with a 
complex amplitude x(t) = -y/8y/(m0 2 )ci'(T)/ y/n, from Eq. (HJi 
we obtain the same EOM for x(r) as above, but with a quan- 
tum mechanically shifted detuning 

A = -2e(<5-3 r /N), (5) 

instead of the classical result A = -2Q6. In classical case, S 
is large, e.g., N 366 as we estimated from the experimental 
circuit parameters, which makes the quantum shift negligibly 
small. Nevertheless, in mesoscopic regime, e.g., S ^ 12 in 
present case, this quantum shift cannot be neglected, as shown 
in Fig. 2(b) and (c), where the critical driving strength moves 
to 0.11F C . Below we illustrate this quantum shift in phase 
space in terms of the Wigner function. 

The Wigner function is defined as: Wix, p, f) — 
l/inh) £™{x + x'\pit)\x - x') e\pi-i2px'/h)dx'. In Fig. 3(a) 
and (b), for comparative purpose, we plot the Wigner function 
simply using p s = \a){a\, where the coherence number a is 
determined from x, based on the steady-state solution of the 
amplitude EOM without and with accounting for the quan- 
tum shift, while in Fig. 3(c) and (d) we show the results from 
direct simulation. We see that in mesoscopic regime it is es- 
sential to account for the quantum shift, in order to make the 
amplitude EOM agree with numerical simulation, as indicated 
by the dashed vertical lines in Fig. 3. (The vertical line in each 
sub-figure indicates the valid "x" center of Wigner function of 
the SAS.) It is then observed that the classical result in Fig. 
3(a) has considerable deviation. Moreover, from Fig. 3(d) we 
get an insight that starting on from certain transient stage, the 
oscillator evolves into a mixed state, which can be formally 
expressed as: Wix,p,t) = P s it)W s ix, p) + P L {t)W L {x, p,t). 
Here, Ws ix, p) and Wl(x, p, t) are the Wigner functions of the 
intrinsic SAS and LAS, say, the two attractors, while Psit) 
and P E U) are the respective occupation probabilities. 

V. TUNNELING RATE AND SCALING EXPONENT 

Based on the structure of the mixed state Wix, p, t), we can 
formulate a way to determine the dynamic tunneling rate from 
SAS to LAS as follows. Originally, in laboratory frame, this 
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FIG. 2: (Color online) (a) and (b): Bifurcation diagram of the DDO 
based on the EOM of x(t), (a) without and (b) with accounting for 
the quantum shift of the effective detuning [see Eq. {5)]. The dashed 
vertical lines indicate the respective bifurcation points, (c): Results 
from numerical simulation of x(t) =Tr[xp(r)] based on Eq. ([3j. We 
choose the x(t) at t = 160 * (2n/Q) to represent the steady state 
amplitude x. The circles and squares stand for results from different 
initial conditions (i.e. the SAS and LAS), from which we observe the 
hysteresis behavior. Here, again, the dashed vertical line indicates 
the bifurcation point in (b), say, 0.77F C , which agrees well with the 
numerical simulation. Parameters: k = 0.01, T=5mK, 6 = 0.065, 
and N = 12. 



rate can be extracted from the occupation probability of the 
SAS via T, = -dP s (t)/dt, while P s (t) = {ae iv <\p(t)\ae iv <). 
Note that here we have used the fact that the SAS in labo- 
ratory frame is a rotating coherent state. More conveniently, 
transformed to the rotating frame, Ps{t) = Tr[p s p(t)], where 
p s = U\f)p s (t)U(t) = \a)(a\ is static, and p(f) = U' (t)p(t)U(t) 
is the DDO state in the rotating frame described by Eq. (J3). 

In the inset of Fig. 4 we show a representative Ps(t), ob- 
tained from numerical simulation using Eq. @ and the for- 
malism outlined above, from which we see that a single expo- 
nential fit can well characterize Ps(t)- This indicates nothing 
but a rate process from SAS to LAS. Moreover, the success of 
a single exponential fit indicates a dominant forward process 
from SAS to LAS. That is, at the early escape stage, the back- 
ward process from LAS to SAS has not yet happened. This 
is similar to the situation in determining the tunneling rate in 
double-well problem, where the backward tunneling process 
is negligibly small at certain early time stage. 

Therefore, using this numerical approach we can extract the 
dynamic tunneling rate (r,) from SAS to LAS. Further, fol- 
lowing Ref. ifioll . we assume V, = Ce" R ^ x . In this proposed 
form, C is an irrelevant prefactor, while the exponential factor 
e -RM originates from an effective activation process. In lim- 
iting cases, such as for classical thermal activation, R is the 
activation energy and A the temperature; while for quantum 
tunneling through a barrier, R is the tunneling action and A the 
Plank constant. In our present case, it is a generalization. We 
may view R as an effective activation energy and A an effective 
Planck constant or temperature. 

Very interestingly, it was found in Ref. ifioll that the dy- 
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FIG. 3: (Color online) (a) and (b): Wigner function of the coherent 
state determined from the steady-state solution of the EOM of x(t), 
(a) without and (b) with accounting for the quantum shift, (c) and (d): 
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namic tunneling action R displays a perfect scaling behavior 
with the driving distance to the critical point, which is defined 
as j] — - F2. Quantitatively, it was found R o: rf and 
a = 3/2. Here, for the mesoscopic DDO, owing to the quan- 
tum shift of the critical point, which moves to 0.77F C under 
the assumed parameter condition as shown in Fig. 2, we define 
accordingly the driving distance as rj — (0.77 F c ) 2 - F?. Re- 
markably, we demonstrate in Fig. 4 that the scaling behavior 
ofRccrj" still holds, yet with an alternative scaling exponent, 
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a 1, instead of a = 3/2 as found by Dykman ifioll . 

We noticed that in Ref. JHfl, scaling behavior of the tran- 
sition rate with the driving frequency (but not the driving 
strength) was analyzed to give a 1.3 ~ 1.4, by a rough 
fitting from a few experimental data. Meanwhile, in the ex- 
periments by Siddiqi et al. i29ll . an effective potential with a 

barrier height scaled as At/° n oc [l - (Fq/F c ) 2 ^ , was em- 
ployed to analyze their measured data. By means of a thermal- 
activation rate oc exp[-AU® vn /k B T^j, this would result in a 
scaling exponent a = 3/2, while a rough quantum WKB es- 
timate for the quantum tunneling rate oc exp^- ^ AU° yn a / hj, 

where a is an effective width of the barrier, is seemingly to 
give or = 3/4. 

Moreover, in Ref. JH], it was shown that both the exponents 
"3/2" and "1" can be obtained, respectively, resulting from 
a local and nonlocal Hamiltonian bifurcation. Nevertheless, 
we should distinguish Ref. Jsij] from our present study in at 
least two aspects: (i) The rates obtained in Ref. JH] are classi- 
cal ones, i.e., the thermal-activation-caused escape rates, since 
they are respectively related to an "effective potential" that is 
then inserted into the thermal-activation rate formula. ( ii) For 
the DDO model under our present study, the "local bifurca- 
tion" is near the larger bifurcation point, and the "nonlocal 
bifurcation" is close to the smaller one. This clarification is 
clearly stated in Sec VII of Ref. JH]. Our numerical study, 
however, is focused on the region near the larger bifurcation 
point, thus should correspond to the "local bifurcation", and 
should be compared with the exponent "3/2". 

In a recent workf2(J, we performed a real time simulation 
for the quantum dynamics of the mesoscopic DDO in labora- 



tory frame, in which the driving field was not taken into ac- 
count in the dissipation terms. Also, in Ref. lEoll . the quantum 
shift was canceled owing to the extra term "x 2 $ / (2m,w 2 )" 
in the Caldeira-Legget model. In the present work, however, 
absorbing this term into the system Hamiltonian, which effec- 
tively renormalizes the DDO's frequency, makes the quantum 
shift evident. In regard to the scaling behavior, quite desirably, 
both studies of Ref. i2(ill and the present one give consistent 
result. This strongly implies that the scaling exponent a = 3/2 
is non-universal, and a — 1 is an alternate scaling exponent in 
the mesoscoic regime. Further investigation for the crossover 
behavior from a mesoscopic to the usual classical regime is 
of particular interest. However, such kind of study will en- 
counter increasing difficulty in real time simulation, which is 
to be the task of our forthcoming research. 



VI. SUMMARY 

To summarize, in a mesoscopic regime we investigated the 
dynamic quantum tunneling of the driven Duffing oscillator. 
Owing to the mesoscopic nature, we found that the critical 
driving strength has a quantum shift, and the tunneling ac- 
tion (extracted exponentially from the tunneling rate) exhibits 
a perfect linear scaling behavior with the driving distance to 
the quantum shifted critical point. 
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